// Gmsh - Copyright (C) 1997-2019 C. Geuzaine, J.-F. Remacle
//
// See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/gmsh/gmsh/issues.
//
// Contributed by Matti Pellikka <matti.pellikka@gmail.com>.

#ifndef HOMOLOGY_H
#define HOMOLOGY_H

#include <sstream>
#include "CellComplex.h"
#include "ChainComplex.h"
#include "Chain.h"
#include "OS.h"
#include "GModel.h"
#include "Options.h"

#if defined(HAVE_KBIPACK)

std::vector<int> vecN0(int n);

// Interface class for homology computation in Gmsh
class Homology {
private:
    // the Gmsh model for homology computation
    GModel *_model;

    // domain and the relative subdomain of the homology computation
    // physical group IDs
    std::vector<int> _domain;
    std::vector<int> _subdomain;
    std::vector<int> _nondomain;
    std::vector<int> _nonsubdomain;
    std::vector<int> _imdomain;
    // corresponding geometrical entities
    std::vector<GEntity *> _domainEntities;
    std::vector<GEntity *> _subdomainEntities;
    std::vector<GEntity *> _nondomainEntities;
    std::vector<GEntity *> _nonsubdomainEntities;
    std::vector<GEntity *> _immuneEntities;

    // save original cell complex
    bool _saveOrig;

    // use cell combining
    int _combine;
    // use cell omit
    bool _omit;
    // use chain smoothning
    bool _smoothen;
    // corecution heuristic
    int _heuristic;

    // file name to store the results
    std::string _fileName;

    // cell complex of the domain
    CellComplex *_cellComplex;

    // whether representatives of (co)homology bases are available
    bool _homologyComputed[4];
    bool _cohomologyComputed[4];

    // resulting betti numbers and chains
    int _betti[4];
    std::vector<Chain<int> *> _chains[4];
    std::vector<Chain<int> *> _cochains[4];

    typedef std::map<Cell *, int, Less_Cell>::iterator citer;

    void _getEntities(const std::vector<int> &physicalGroups,
                      std::vector<GEntity *> &entities);
    void _getElements(const std::vector<GEntity *> &entities,
                      std::vector<MElement *> &elements);

    // create a string describing the generator
    std::string _getDomainString(const std::vector<int> &domain,
                                 const std::vector<int> &subdomain) const;

    // construct the cell complex
    void _createCellComplex();

    // create and store a chain from homology solver result
    void _createChain(std::map<Cell *, int, Less_Cell> &preChain,
                      const std::string &name, bool co);

    void _deleteChains(std::vector<int> dim = vecN0(4));
    void _deleteCochains(std::vector<int> dim = vecN0(4));

    std::vector<int> _addToModel(int dim, bool co, bool post,
                                 int physicalNumRequest) const;

public:
    // Determine domain and relative subdomain of (co)homology computation
    Homology(GModel *model, const std::vector<int> &physicalDomain,
             const std::vector<int> &physicalSubdomain,
             const std::vector<int> &physicalIm, bool saveOrig = true,
             int combine = 3, bool omit = true, bool smoothen = true,
             int heuristic = 1);
    ~Homology();

    GModel *getModel() const { return _model; }
    void setFileName(const std::string &fileName) { _fileName = fileName; }

    void getDomain(std::vector<int> &domain) const { domain = _domain; }
    void getSubdomain(std::vector<int> &subdomain) const
    {
        subdomain = _subdomain;
    }

    // find the bases of (co)homology spaces
    // if dim is not provided,, find 0-,1-,2-,3-(co)homology spaces bases
    // otherwise only find those indicated in dim
    void findHomologyBasis(std::vector<int> dim = vecN0(4));
    void findCohomologyBasis(std::vector<int> dim = vecN0(4));

    // find a homology and cohomology basis pair such that
    // the incidence matrix of the bases is an identity matrix
    // if master==0, homology basis determines the cohomology basis
    // if dim is not provided, find 0-,1-,2-,3-(co)homology spaces bases
    // otherwise only find those indicated in dim
    void findCompatibleBasisPair(int master = 0, std::vector<int> dim = vecN0(4));

    // is the (co)homology in given dimensions already compited
    // if dim is not provided, return true only if computed in all dimensions
    bool isHomologyComputed(std::vector<int> dim = vecN0(4)) const;
    bool isCohomologyComputed(std::vector<int> dim = vecN0(4)) const;

    // add chains to Gmsh model
    // dim: only add dim-chains if dim != -1
    // post: create a post-processing view
    // physicalNumRequest: number the chains starting from this if available
    // returns physical group numbers
    std::vector<int> addChainsToModel(int dim = -1, bool post = true,
                                      int physicalNumRequest = -1) const;
    std::vector<int> addCochainsToModel(int dim = -1, bool post = true,
                                        int physicalNumRequest = -1) const;

    // get representative chains of a (co)homology space basis
    void getHomologyBasis(int dim, std::vector<Chain<int> > &hom);
    void getCohomologyBasis(int dim, std::vector<Chain<int> > &coh);

    // find Betti numbers
    // faster than finding bases for (co)homology spaces
    void findBettiNumbers();

    // are Betti numbers computed
    bool isBettiComputed() const;

    // get a Betti number
    int betti(int dim);

    // get the Euler characteristic
    int eulerCharacteristic();

    // write the generators to a file
    bool writeBasisMSH(bool binary = false);

    // store dim-dimensional cells of cellComplex as a physical group
    // in _model, for debugging
    void storeCells(CellComplex *cellComplex, int dim);
};

#endif

#endif
